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Abstract 

We use cosmological gas dynamic simulations to investigate the accuracy of galaxy cluster 
mass estimates based on X-ray observations. The experiments follow the formation of clus- 
ters in different cosmological models and include the effects of gravity, pressure gradients, 
and hydrodynamical shocks. A subset of our ensemble also allows for feedback of mass 
and energy from galactic winds into the intracluster medium. We find that mass estimates 
based on the hydrostatic, isothermal /3-model are remarkably accurate when evaluated at 
radii where the cluster mean density is between 500-2500 times the critical density. Applied 
to 174 artificial ROSAT images constructed from the simulations, the distribution of the 
estimated-to-true mass ratio is nearly unbiased and has a standard deviation of 14-29%. 
The scatter can be considerably reduced (to 8-15%) by using an alternative mass estimator 
that exploits the tightness of the mass-temperature relation found in the simulations. The 
improvement over /3-model estimates is due to the elimination of the variance contributed 
by the gas outer slope parameter. We discuss these findings and their implications for 
recent measurements of cluster baryon fractions. 
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1. Introduction 



Estimates of the total masses of groups and clusters of galaxies have been used to infer the 
amount of dark matter in the universe for over sixty years. A novel technique for estimating 
the density parameter fio makes use of precise measurements of the visible, baryonic mass 
fraction /& in galaxy clusters along with limits on the universal baryon fraction derived 
from primordial nucleosynthesis. White et al. (1993) showed that cluster baryon fractions, 
defined as the ratio of the mass in galaxies and intracluster gas to the total cluster mass, 
should not differ substantially from the universal value, VL^/VLq, when determined near the 
outer boundary of their hydrostatic regions — roughly an Abell radius for a cluster as rich 
as Coma. A straightforward and unbiased estimate of the density parameter is formed by 
Q = Q b /f b . 

Recent measurements in rich clusters (Briel, Henry & Bohringer 1992; Durret et al. 1994; 
David, Jones & Forman 1995; White & Fabian 1995) and poor clusters and groups (Pon- 
man et al. 1994; Pildis, Bregman & Evrard 1995; Neumann & Bohringer 1995) indicate 
fb > 0.04 /i~ 3 / 2 . (Hereafter we write Hubble's constant as 100 h km s -1 Mpc -1 .) Taken 
with the nucleosynthesis determination ~0.0125/i -2 (Walker et al. 1991), this implies a 
rather small value of the density parameter, Oq < 0.3/i -1 / 2 , unless primordial nucleosyn- 
thesis calculations have underestimated the universal baryon fraction by almost a factor 
of three. Although there is current debate about the uncertainty in primordial nucleosyn- 
thesis determinations of O fe (Krauss & Kernan 1994; Copi, Schramm & Turner 1995; Hata 
et al. 1995; Steigman 1995; Sasselov & Goldwirth 1995), current interpretation of the data 
appear to rule out the large values of required for consistency with a universe with 
closure density. 

These upper limits on fi are especially strong because most of the baryons in clusters are 
in the hot intracluster medium (ICM), a component empirically found to be more extended 
than the dark matter distribution (e.g., David et al. 1995). Numerical simulations show 
that this is a general result of hierarchical scenarios, where clusters are formed through 
mergers of protoclusters during which energy is transferred systematically to the gas from 
dark matter (Navarro & White 1993; Pearce, Thomas & Couchman 1994). The results 
imply that cluster baryon fractions, measured at radii encompassing a density contrast of 
a few hundred, should be about 10 — 20% lower than the universal value (Evrard 1990; 
Cen & Ostriker 1992). 

Several solutions have been proposed to rescue = 1 from what has been deemed the 
cluster "baryon catastrophe" (Carr 1993). A low Hubble constant H < 30 km s _1 Mpc -1 
can alleviate the problem (Bartlett et al. 1995), but is in strong disagreement with recent 
observational estimates which favour a high value for H (Freedman et al. 1994; Mould 
et al. 1995). A simple possibility which has not yet been explored in detail is that the 
binding masses inferred from X-ray observations may be systematically underestimated 
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by a significant amount. These estimates usually rely on assumptions such as spherical 
symmetry, hydrostatic equilibrium for the gas, and virial equilibrium for the galaxies, 
none of which may be fully realized in practice. Clusters exhibit signatures of substructure 
both in the ICM (Mohr et al. 1995; Buote & Tsai 1995) and in their galaxy distributions 
(Dressier & Schectman 1988; Bird 1995; Crone, Evrard & Richstone 1995), suggestive of 
recent mass accretion and of significant departures from equilibrium. A moderate bias in 
the mass estimates introduced by these effects would have important consequences for the 
Q limits derived by the above arguments. 

Discrepancies of ~ 50% have, in fact, been reported when mass estimates of the central 
regions of clusters derived from X-ray observations are compared with those required to 
produce strong arcs by gravitational lensing of background galaxies (Miralde-Escude & 
Babul 1995). Although this discrepancy could signal significant departures from hydro- 
static equilibrium or support from non-thermal sources such as magnetic fields (Loeb & 
Mao 1994), systematic errors in the lensing mass estimates due to projection effects and 
substructure are more likely to be responsible for the disagreement (Bartelmann, Steinmetz 
& Weiss 1995; Bartelmann 1995). In fact, at larger radii, weak gravitational lensing has 
also been used to measure cluster masses (Tyson, Valdes & Wenk 1990; Bonnet, Mellier 
& Fort 1994; Fahlman et al. 1994) and the small number of clusters with both X-ray and 
weak lensing mass estimates show no significant discrepancy between the two (Smail et al. 
1995; Squires et al. 1995), although the statistical uncertainties are still large. 

In this paper, we examine the accuracy of X-ray binding mass estimates using a large 
number of high resolution simulations designed to follow the non-linear, dynamical evo- 
lution of the gravitationally coupled system of dark matter, gas and (in one set of runs) 
galaxies in a variety of cosmological models. A total of 58 clusters are drawn from three 
separate projects using two different Lagrangian hydrodynamical codes. The models are 
used to generate 174 synthetic ROSAT X-ray images and broad beam temperature esti- 
mates. Binding masses are then estimated for these systems in a manner analogous to 
that applied to observational datasets. These models are therefore ideal for understanding 
possible systematic effects afflicting cluster mass estimates based on X-ray observations. A 
recent paper by Schindler (1995) addresses the same issue using a sample of six simulations 
generated with very different techniques. Our results are in excellent agreement, despite 
the fact that the numerical methods and data analysis procedures used in the two studies 
differ in a number of details. 

After describing briefly the numerical simulations in §2, we begin by examining the validity 
of the hydrostatic and isothermal assumptions using the three dimensional velocity and 
temperature profiles of simulated clusters (§3). We then investigate the accuracy of bind- 
ing masses estimated using the simplest combination of X-ray imaging and broad beam 
temperatures (§4). In §5 we discuss how the tight correlation between cluster mass and 
X-ray temperature can be used to determine binding masses with even greater statistical 
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accuracy. We conclude in §6 with a brief discussion of some implications of these results. 

2. Numerical Methods 

2.1 Sample Description 

We use 58 N-body/gas dynamics simulations drawn from three different projects and run 
with two completely independent Lagrangian codes. In all cases we use the Smoothed 
Particle Hydrodynamics (SPH) technique to follow the evolution of the gas, and either a 
P 3 M code or a tree-based N-body code to compute the gravitational interaction between 
particles. All the simulations assume a standard, cold dark matter (CDM) initial fluctu- 
ation spectrum with F = flh — 0.5. We neglect the radiative cooling of the gas, as well as 
magnetic fields as a possible source of pressure support. We take the baryon density pa- 
rameter to be 0.1, and use a Hubble constant of 50 km s -1 Mpc -1 (h = 0.5) when scaling 
to physical units. Details of the models can be found in the references quoted below. 

The first set consists of 28 runs from Chris Metzler's thesis (Metzler 1995), which examines 
the effects of energy feedback and mass ejecta from early type galaxies on the evolution of 
the intracluster medium (ICM). The set consists of 14 clusters of different mass obtained 
by constraining the initial density field in cubic, periodic regions ranging from 25 to 60 
Mpc (Bertschinger 1987). Each realization is run twice, once including the effects of galaxy 
feedback and a second, control run without feedback. We shall refer to each of these samples 
as "EJ" and "2F", respectively. All of the runs assume = 1. Each simulation uses a 
total of 65, 536 particles, divided equally between gas and dark matter. A description of 
the feedback implementation and application to the formation of a Coma-sized cluster can 
be found in Metzler & Evrard (1994). The full set of runs is described in Metzler & Evrard 
(1995). A salient feature of the ejection runs is that they employ a rather extreme model of 
galactic feedback, in which early-type galaxies lose half their initial mass by winds. The EJ 
and 2F series are intended to define an envelope within which realistic models of feedback 
should lie. 

The second set of runs is a sample of 24 used to investigate the X-ray cluster morphology- 
cosmology connection by Evrard et al. (1993) and Mohr et al. (1995). For this project, 
eight different realizations of the initial density field were evolved within three background 
cosmologies in periodic cubes ranging from 30 to 60 Mpc in length. Again, 65, 536 particles 
per run were used. The cosmological models explored were the standard CDM scenario 
(model "EdS"; as = 0.59, O = 1), an unbiased, open CDM universe (model "Op2"; 08 = 1.0, 
Oo = 0.2, A o = 0) and an unbiased, low density universe with a flat geometry (model "F12"; 
cr 8 = 1.0, fi = 0.2, A = 0.8). Here cr 8 is the rms mass fluctuation in spheres of 8/i _1 
Mpc, and O and Ao are the present values of the density parameter and the cosmological 
constant, respectively. An important aspect of this set is that the runs corresponding 
to each cosmology share common dynamical histories through the use of the same eight 
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initial density fields. The benefit of this procedure is that relative differences between the 
final characteristics of clusters can be ascribed to the effect of the different cosmological 
backgrounds rather than to "cosmic scatter" in the initial density fields. 

The final set consists of six models taken from Navarro, Frenk & White (1995a; model 
"NFW" ) which were used to examine the evolution of scaling laws relating the dynamical 
and X-ray properties of clusters of different mass. These simulations were evolved with a 
code completely independent from the one used in the runs described above; a Tree/SPH 
code described in detail in Navarro & White (1993, 1994). The six clusters were identified 
at z = (a 8 = 0.63) in the cosmological N-body simulations of Frenk et al. (1990), and 
then resimulated individually at higher resolution. Each simulation has 21,296 particles, 
half gas and half dark matter, distributed over boxes of size 15-50 Mpc, depending on the 
mass of the cluster. The tidal field due to material surrounding each cluster out to 360 
Mpc is treated self-consistently by coarse sampling the surrounding matter with particles 
of radially increasing mass. 

The ensemble of simulations span a wide range in mass and temperature, from ~ 10 14 to 
3 x 1O 15 M , and from ~ 1 to 10 keV. Table 1 provides a summary for future reference of 
the runs and the notation that we use. The clusters produced in the different projects have 
similar spatial and mass resolution properties. The ratio between the size of a simulated 
cluster (as measured by the radius, rsoo, where the mean density relative to the critical 
value is 500) and the gravitational softening is in the range 10 < r$oo/e < 30. The runs 
have similar numbers of particles within r5oo, typically 5000 in each component. 

2.2 Hydrostatic, Isothermal Mass Estimates 

The typically smooth morphology of the X-ray emission from the hot, intracluster medium 
leads naturally to the hypothesis that the gas is near equilibrium, stratified along isopo- 
tential surfaces in a mildly evolving distribution of dark matter, gas and galaxies. The as- 
sumption of hydrostatic equilibrium — the balance between pressure gradients and gravity 
- for gas supported solely by thermal pressure, results in a direct measure of the binding 
mass M(r). Assuming spherical symmetry, 

M(r) = kT ( r ) r ( dl °Sp(r) | d log T(r) \ 
Gfirrip \ d log r d log r / 

where p(r) and T(r) are the gas density and temperature profiles, k is Boltzmann's con- 
stant, and \xm v is the mean molecular weight of the gas. In principle, all terms on the 
right hand side of this equation are measurable. The main limitation is that one must 
deconvolve three dimensional profiles from two dimensional surface brightness informa- 
tion. This requires knowledge of, or a model for, the temperature profile T(r). Since 
direct measurements of the temperature as a function of radius T(r) (more precisely, the 
X-ray emission-weighted projected temperature) are still a relatively rare commodity, the 
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common practice has been to assume that the gas is isothermal at the spatially averaged 
temperature Tx determined from a broad beam spectroscopic instrument such as EXOSAT 
or Ginga. Preliminary ASCA results (Ikebe et al 1994; Mushotzky 1994), as well as pre- 
vious direct measurements (see, e.g., Watt et al. 1992) and the numerical models that we 
use here generally support this assumption (see §3.3 below). 

The usual parametrization of the density profile of the ICM is based on the isothermal, 
/3-model proposed by Cavaliere & Fusco-Femiano (1976), Picm (r) = / 9 (l + (r/r c ) 2 )- 3 ^/ 2 . 
With this assumption, equation (1) reduces to 

Mir)- W kTxr (r/rc)2 - 1 13 x 1Q 15 3 Tx r (r/rc)2 M, , (2) 
M[r) - G ^m p l + (r/r c )2 - L13X1 ° ? 10 keV Mpc 1 + (r/r c ) 2 Mq ' (2) 

assuming \i = 0.59 for the second equality. For an isothermal gas, the X-ray surface 
brightness as a function of projected radius Sx(0) simply reflects the integrated emission 
measure and can be expressed as 

S X (9) = So (l + (6/6 c ) 2 )- 3P+1/2 - (3) 

X-ray spectroscopy provides Tx, while estimates of 3 and r c are obtained from X-ray 
imaging. 

To mimic the mass estimates derived from X-ray observations, we generate for each simu- 
lated cluster three artificial ROSAT images and emission weighted temperatures Tx along 
the principal axes of the volume. Generating synthetic ROSAT images requires a choice 
of several parameters, including: (i) the cluster redshift z, (ii) the exposure time t exp , and 
(iii) the background noise level Sn- Fitting the resulting surface brightness profiles to 
the /3-model introduces two additional parameters; (iv) the background subtraction level 
Sb and (v) the minimum flux level to which the fit is performed S m i n . An additional 
choice is the energy band of the observations, which we take to be the full ROSAT band 
0.1-2.4 keV. Because there is no obscuring galaxy in our 'observations', our results are 
insensitive to the exact choice of the lower energy limit for the band employed. We use 
emission weighted temperatures within circular apertures large enough to contain nearly 
all the cluster emission. Since most of the photons come from the central regions of the 
cluster, the estimates of Tx are rather insensitive to the aperture choice. 

We use a single parameter set to construct the synthetic images of all model clusters, 
chosen to be representative of typical X-ray observations. We view the clusters at a 
redshift z = 0.04, image them for t exp = 7200 s, include a Poisson background in the 
counts at a level of 3 x 10~ 4 cnts s -1 arcmin -2 , subtract a mean background with exactly 
this value, and fit each to the /3-model out to the radius at which the corrected counts 
reach the background level. We use emission weighted temperatures in the same ROSAT 
band, measured within an angular scale of 16', corresponding to a fixed metric radius of 
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1.04 Mpc. These parameter values are summarized in Table 2. We investigated alternative 
parameter sets, and found that the results depended mainly on the image quality through 
the combination of redshift, exposure time and assumed background. The parameter 
choices above provide high signal-to-noise images for the majority of the sample. 

For each cluster, the three images and emission weighted temperature maps are generated 
and reduced to obtain the values of /3, r c and Tx. This procedure works very well for 
clusters formed in cosmological models with = 1. However, as discussed in Mohr et al. 
(1995), in low density universes the emission profiles of the cluster models are strongly 
peaked, resulting in unacceptably poor /3-model fits. In this case, we decided to excise the 
central regions and to fit only data outside a region of angular radius 4', corresponding to 
a linear scale of 260 kpc (see Mohr et al. 1995 for details). Typical statistical uncertainties 
in the fitted parameters are < 5%, but in a few cases the uncertainties can be as high as 
~40%. Once r c and T x have been determined, we can estimate the binding mass as a 
function of radius using equation (2). 

2.3 A Particular Example 

Before examining results for the entire sample, we discuss data for a particular cluster 
which highlight characteristics typical of the ensemble. Figure 1 shows X-ray surface 
brightness and emission weighted temperature maps, as well as the emission weighted 
projected velocity field for the ICM in a cluster taken from the EdS sample. The field of 
view in each of the panels is 64', corresponding to 4.2 Mpc at the assumed cluster redshift 
of 0.04. Two of the three projections show a bimodal central structure; the result of a recent 
merger involving two main sublumps with mass ratio 2.8:1. Infall patterns of the sublumps 
are evident in the velocity field of the middle (y-axis) projection. The temperature map 
shows significant variations (up to a factor ~ 2) in temperature on scales of a few hundred 
kpc. The hot spots occur in regions where the gas is being compressed and mildly shocked 
by the interpenetrating subcluster cores. Cooler gas can be seen trailing in the wake of 
the substructure cores. 

Several of the features of this simulated cluster, particularly the y-axis view, are reminiscent 
of the cluster A2256. Features in common include a bimodal central structure and spatial 
variations in temperature similar in morphology and amplitude (Briel et al. 1991; Miyaji 
et al. 1993; Briel & Henry 1994). Furthermore, an offset, extended radio halo exists in 
A2256 which strongly indicates the presence of a flow pattern similar to that seen in the 
y-axis image (Rottgering et al. 1994). 

Figure 2 shows the azimuthally averaged X-ray surface brightness and emission weighted 
temperature profiles for the three projections. Surface brightness profiles are scaled down 
successively by an order of magnitude for the sake of clarity; as are the temperature profiles 
by factors of 2. Values of the fitted parameters are r c = 294 ± 6, 384 ± 10, and 359 ± 9 
kpc; and /3 = 0.89 ± 0.02, 0.85 ± 0.02, and 0.81 ± 0.02 for the three projections from top to 
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bottom, respectively. 

Although significant spatial variations in temperature are obvious in Figure 1, the radially 
averaged temperature profile is close to isothermal over a significant fraction of the cluster 
image. The temperature varies by < 20% within 10', a region wherein the surface bright- 
ness drops by a factor of 30. The temperature drops by a factor of 2 at about 20' from the 
center, but the surface brightness at this radius is already below the adopted background 
of 3 x 10 -4 cnts s _1 arcmin -2 . 

From Figure 1 we learn the importance of examining spatial temperature maps directly, 
since a flat azimuthally averaged profile need not be indicative of a truly isothermal ICM. 
Comparison of the temperature and surface brightness maps can provide useful dynamical 
clues, although geometry plays an obscuring role. From the single dynamical configuration 
corresponding to the cluster in Figure 1 one can get, depending on projection (from top to 
bottom), (i) a fairly relaxed X-ray image with asymmetric temperature map and strong 
velocity gradients; (ii) a bimodal X-ray image with an "S" -shaped hot spot resulting from 
a symmetric infall pattern and (iii) a bimodal X-ray image with a peanut-shaped hot spot 
and a relatively modest velocity field in projection. 

3. Is the ICM Hydrostatic and Isothermal ? 

The assumption of hydrostatic equilibrium underlies the X-ray mass estimate method, and 
the assumption of isothermality greatly simplifies it. In this section, we analyze how close 
to hydrostatic equilibrium and isothermality the cluster models are by inspecting directly 
their three dimensional velocity and temperature fields. 

3.1 The ICM Velocity Field 

Figure 3 shows the radial Mach numbers (v r )/c s and (v 2 )/c 2 s derived from the gas sound 
speed c 2 s {r) =5kT(r)/3/j,m p and the first and second moments of the radial velocity field, 
respectively. To facilitate comparison between clusters of different sizes, the mean interior 
density contrast 8 c (r) =3M(r)/47rp c r 3 is used as the radial variable in the figure, reversed 
to reflect the correspondence with cluster radius. Note that the density contrast used here 
is defined with respect to the critical value for closure p c = 3 Ho 2 /87rG in all the models. 
The center of the cluster is defined as the position of the most bound dark matter particle, 
and velocities are calculated with respect to the mean velocity of all cluster particles linked 
by a standard friends-of-friends algorithm using a linking parameter 0.15 times the mean 
interparticle spacing. The = 1 runs all exhibit a common structure, with the exception 
of two of the NFW runs and one EdS run, which are undergoing strong mergers at z = 0. 
The Mach numbers of these systems are significantly higher than the average. The low O 
runs have generally quieter velocity fields, as expected given their earlier formation times 
and, consequently, their dynamical maturity relative to clusters formed in a high density 
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universe. 



The mean radial Mach number (v r )/c s has the characteristic signatures expected from 
gravitational infall (Gunn & Gott 1972). An outer zone of mildly supersonic infall sur- 
rounds the 'virial' region of the cluster, within which the gas has been largely thermalized 
and is close to hydrostatic equilibrium. The infall regime is largely absent in the low O 
models, due to the stagnated growth of linear perturbations on large scales. This, however, 
does not imply that there are no recent merger events in the open universe sample — at 
least one Op2 cluster is experiencing an ongoing merger event at the present time. 

The right-hand panels of Figure 3 provide an upper limit on the ratio of kinetic to thermal 
pressures for the gas. This ratio rises monotonically with radius, from values < 10% for 
radii where 8 C > 500 to values > 50% at radii where S c < 100. A few ongoing mergers are 
clearly recognized by the large values of the {v^)/c 2 s ratio near the center. In this case, 
interpreting the ratio between kinetic and thermal energy as the relative contribution of 
kinetic and thermal pressures to support the gas does not apply, since the systems are far 
from equilibrium. As the velocity field in Figure 1 indicates, non-zero values of (v%) arise 
during mergers from large-scale bulk motions of the gas across the face of a given radial 
shell rather than from a local, uniform dispersion on the shell. 

From Figure 3, we derive a value S c = 500 as a conservative estimate of the boundary 
between the inner, virialized region of the clusters and their recently accreted, still settling 
outer envelopes. Define r§ c as the radius within which the mean interior density is 5 C times 
the critical value. Then, within rsoo, the hydrostatic equilibrium assumption is valid since 
the gas is, on average, neither expanding nor contracting. This estimate is conservative in 
the sense that, in many clusters, hydrostatic balance appears to hold even at somewhat 
larger radii. The turnaround radius can be seen to occur at radii about a factor 3 — 4 larger 
than r 5 oo (for = 1). This is consistent with the spherical infall models of Bertschinger 
(1985) in which infalling gas is shocked nearly to rest at a radius about a third of the 
turnaround radius. Despite this nice agreement, we stress that the accretion in these three 
dimensional models is far from spherically symmetric. 

Table 3 gives mass averaged values of the two radial Mach numbers measured within r 50 o 
for each set of runs. The mean radial Mach numbers are all quite small, typically a few 
percent or less, and are consistent with zero given the measured error in the mean. Again, 
the NFW set seems to be the most dynamically active, as can be seen from the measures 
of both velocity moments. Typical values of {vf)/c 2 are < 10%, indicating that the gas is 
hydrostatically supported by thermal pressure to this accuracy within rsoo- 
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3.2 Cluster Scaling 



As discussed above, we will use the radius rsoo as a characteristic length scale separating the 
nearly hydrostatic central region of a cluster from the surrounding, recently accreted outer 
envelope. As shown by Navarro et al. (1995a,b) and Metzler & Evrard (1995), clusters of 
different mass have similar structures when scaled to such a characteristic radius. This 
similarity, along with the equilibrium assumption validated above, implies a power-law 
relationship between cluster 'size' and temperature 

M(< r 50 o) 2 
T oc oc r 500 , (4) 

?~500 

as appropriate for systems of similar density in virial equilibrium. Figure 4 shows this 
relationship for the simulations, using the measured value of rsoo an d three orthogonal 
measures of the emission weighted temperature T x . Each cluster appears three times in 
this plot; the dispersion in Tx for different projections is typically quite small. Fitting the 
results with a power law, we find 

r 50 o(T x ) = 2.48 ±0.17 (— ^) ' Mpc (5) 

where the quoted uncertainty in the intercept is given by the standard deviation of the 
residuals in the log space fit. The actual best fit slope for each individual set of runs differs 
by less than 10% from the 0.5 exponent expected from equation (4). There are small 
offsets between the models. The ejection models, for example, have ~6% smaller rsoo at 
a given temperature or, equivalently, ~ 12% higher emission weighted temperatures for a 
given rsoo, compared to the 2F runs. The slope for the EJ sample is also slightly (~10%) 
steeper than those of the non-ejection samples, as expected from the differential effect of 
feedback, which raises the temperature of poor cluster gas proportionally more than that 
of rich clusters (Metzler & Evrard 1995). There is no significant difference between the 
NFW, 2F and EdS data sets, indicating an encouraging agreement between the results of 
models run with completely independent codes for the same cosmological model. 

3.3 Temperature Profiles 

Figure 5 shows the mass- weighted gas temperature, in units of Tx, as a function of the 
normalized radius r/rsoo for all runs. For = 1, the profiles are close to isothermal within 
^50o/2, and decline gently beyond that radius; the temperature at rsoo is ~20% lower than 
at the center. In all cases, the modest drop in temperature within rsoo is due to the fact 
that the density profiles of both gas and dark matter are slightly steeper than isothermal 
in their outer parts (Navarro et al. 1995a,b). In the case of an open universe, the profiles 
are noticeably steeper; the temperature at rsoo is on average a factor of two lower than 
the central value. This comes as no surprise, for the density profiles of clusters formed 
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in low-density universes are expected to be significantly steeper than those formed in an 
= 1 cosmology (Hoffman 1988, Crone, Evrard & Richstone 1994). Since X-ray data 
rarely extend significantly beyond rsoo, an isothermal assumption should be appropriate 
for most observed clusters. 

4. Mass Estimates from the /3-model. 

After constructing the synthetic ROSAT images as described in §2.2, fitting the surface 
brightness with equation (3), and measuring the projected X-ray emission weighted temper- 
ature, we derive estimates of the cluster binding mass as a function of radius via equation 
(2). Once again, in order to compare clusters of different mass (and size), we transform 
the radial variable into an estimated density contrast 6^, st relative to the critical density 
using the estimated mass M est (r) from eq.(2); 

zest M = 3M"«(r) 2GM^(r) (2a(r)\ 2 
c K } ~ Anp c rZ H^r 3 ~ V H r J W 

where cr 2 (r) = GM est (r)/2r is roughly the one-dimensional velocity dispersion of the 
cluster. For a rich cluster with a = 1000 km s _1 at r = 1 h~ x Mpc, the estimated density 
contrast is 5^* = 400. Note that b e c st depends on the combination H r, and therefore it is 
independent of the Hubble constant, making it a useful measure of radius for observations. 

Figure 6 shows the ratio between estimated and true mass as a function of the radial 
coordinate 5^ st . The finite dynamic range in the simulations limits the density contrast 
to values < 5000. This figure shows that the cluster binding masses are on average quite 
accurately determined at overdensities between 500 and 2500. In the outer regions, where 
51 st < 200, masses are typically overestimated because the estimated mass, assumed to 
increase linearly with radius, increases with radius faster than the true mass in this region. 
The effect is more pronounced in the low density runs, where the cluster density profiles are 
steepest. Overestimates by factors up to 3 are seen in the low-density runs at 5^ st = 100. 

Figure 7 presents histograms of the estimated-to-true mass ratio at radii corresponding to 
fiest _ 2500, 500 and 100. These three values of the density contrast sample dynamically 
different regions and span a range comparable to that of current X-ray observations. 
Dashed vertical lines in the figure show ±40% error for reference. The trend toward 
overestimates at low density contrasts noted in Figure 6 is apparent in the rightmost 
column. The omission of the NFW runs at this contrast is technical in origin; these 
simulations do not extend reliably to these low overdensities. 

Relatively large (factor ~ 2) underestimates occur in a few of the = 1 clusters. Three 
of the worst offenders arise from images of a single, strongly bimodal EdS cluster which 
is currently experiencing a major merger. Indeed, ongoing mergers are responsible for 
six of the worst underestimates at 5^ st = 500. Synthetic ROSAT images of these cases 
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are shown in Figure 8. Bold and light circles indicate the estimated and true values of 
rsoo- The complex, multi-peaked structure of the X-ray emission in these images is a 
strong signal of dynamical unrest. Although such objects would probably be rejected by 
observers attempting to apply an equation based on hydrostatic equilibrium, we include 
these cases in our analysis below since they only represent a small fraction of the total 
population. 

The similarity between the = 1 sets suggests that co-adding the runs is appropriate 
in order to compute the ensemble statistics. Figure 9 shows the histogram of estimated- 
to-true mass ratios, X = M est /M true , evaluated at 8^ st = 500 for the 126 images of 
the combined EJ, 2F, NFW and EdS sets. The histogram is nearly Gaussian with mean 
A? = 1.02 and standard deviation ax = 0.29. At this density contrast, the /3-model estimates 
are unbiased and have rather modest scatter. We repeated this procedure at contrasts 
51 st = 2500, 1000, 250 and 100. The means and standard deviations of the resulting 
estimated-to-true mass ratio distributions are given in Table 4 for the = 1 and Oo = 0.2 
ensembles. Values in this Table reflect the trends apparent in Figures 6 and 7; the bias 
and the variance in the mass estimator both increase with increasing radius (i.e., towards 
lower density contrasts). The increase in the variance with radius is probably linked to 
the longer dynamical times of the outer regions. 

Although the estimator does worse in the case of strongly bimodal clusters (see Figure 8) , 
there are cases where the geometry of the projection, together with the interplay between 
the estimated values of /?, r c , and Tx, can result in accurate mass estimates even for 
clusters with suspicious looking X-ray images. As an example, consider the cluster shown 
in Figure 1. The 8^ st = 500 mass estimates for the three orthogonal projections yield 
(from top to bottom) # = 1.17, 1.02 and 0.97. The most symmetric X-ray image incurs 
the largest error while the two images displaying core bimodality are more accurately 
determined. Note that rsoo for this cluster is 1.62 Mpc, or 24' in the figure, well beyond 
the core region. The larger value in the top projection compared to the others is due 
to slightly larger values of (3 and Tx and a smaller core radius compared to the other 
projections. These values result from the fact that the line of sight in the top panel is 
nearly parallel to the collision axis of the penetrating cores. 

In summary, we find /3-model mass estimates to be nearly unbiased and accurate to a few 
tens of percent in the regime 250 < 5^. st < 2500 for the = 1 models and 5^, st > 1000 for the 
Oo = 0.2 sample. A bias toward overestimating masses exists at low values of S^, st . Clusters 
with strongly bimodal or more complex images involve the largest mass underestimates. 
Because of the interplay of T x , (3 and r c , there is not a simple, general connection between 
the properties of the X-ray image and the accuracy of the mass estimates obtained with 
the /3-model. 
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5. Estimates Based on Cluster Scaling Relations. 



As discussed in §3.2, massive clusters within a given cosmology exhibit a remarkably similar 
structure when scaled to a fixed density contrast. Together with the condition of virial 
or hydrostatic equilibrium, this implies that the temperature of a cluster is, on its own, 
a good indicator of the size and mass of the system. Figure 4 shows this result very 
clearly. The tight correlation shown in this figure between Tx and rsoo implies a similarly 
tight correlation between mass and temperature since, by definition, the mean density 
within rsoo is 500 times the critical density. Armed only with the broad beam temperature 
measure Tx, we can thus form an estimate 

3 /2 

M 5 ^(T X ) = 500^p c (r 50 o) 3 = 2.22 x 10 15 M Q (7) 

for the mass within the radius r^oo(Tx) given by equation (5). 

We compare this mass estimate with the true mass within r 500 (Tx) for the = 1 ensemble 
in Figure 10. (Results for the O = 0.2 sets are similar.) The distribution of estimated-to- 
true mass ratios is nearly Gaussian with a standard deviation of only 15%. Note that, in 
this case, no cluster in the = 1 sets has its mass over or underestimated by more than 
40%, regardless of its dynamical state. This procedure can be extended to other values of 
the density contrast in a straightforward way. For a given 5 C , we compute the characteristic 
radii rs c from the numerical sample and fit them to a relation of the form 

1/2 

r«.PW = r 10 (« (8) 

where the normalization ri (5 c ) is the average radial scale of 10 keV clusters at density 
contrast S c . 

The resulting distributions of X are unbiased by construction and have standard deviation 
ax- Table 5 shows the characteristic radii rio(d c ) and the scatter in the mass estimator ax 
for density contrasts 5 C =100, 250, 500, 1000 and 2500. Slight offsets in the characteristic 
radii are evident between the high and low O cosmologies, consistent with the difference 
in cluster density profiles. As in the /9-model estimates, the scatter in the mass estimates 
increases with radius. At density contrasts of a few thousand, the uncertainty in the mass 
estimates is extremely small ax ^ 10%. 

Another way to interpret these results is to consider that the scaling law mass estimate 
is consistent with the hydrostatic, /3-model estimate, equation (2), when evaluated at 
^500 (Tx.) with a characteristic value of (3 given by 

0* = 0.79 (1 + (r c /r) 2 ). (9) 
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The second term is a typically small (< 5%) correction at 5 C = 500. This value for (3 
compares well with measured values of (3 for many well studied, rich X-ray clusters, among 
them Coma (Hughes 1989) and A2256 (Briel et al. 1992). 

Why are the scaling law estimates more accurate than those of the /9-model? Consider the 
variance in the /3-model mass estimated at a fixed radius 

In a universe filled with perfectly hydrostatic, self-similar clusters, all clusters would follow 
a r<5 c (Tx) relation like equation (8) exactly and all would have a fixed value /?* for their 
outer profile slope. Introduction of perturbations in temperature and density off this 
perfect sequence will lead to a non-zero variance in the mass estimate. The only way to 
retain perfect mass estimates is to introduce correlated perturbations A/3//3 = — AT X /T X 
in density and temperature, so that hydrostatic equilibrium is maintained. Uncorrelated 
perturbations in (3 and Tx will lead to a larger variance than that arising from perturbations 
in one parameter alone. 

Figure 11 shows the perturbations measured directly in the simulations. Perturbations in 
(3 are defined with respect to the average, /?* = 0.79, while those in temperature are defined 
with respect to the mean radius-temperature relation, using the known value of rsoo to 
define the unperturbed cluster temperature for each cluster. It is clear from this Figure 
that the data exhibit no correlation between A/3//3 and AT X /T X . The larger variance in 
the /3-model compared to the scaling law mass estimates can thus be understood as arising 
from an additional, independent source of error in the /3-model estimator; the introduction 
of {3 from X-ray imaging is essentially adding noise to the mass estimates. As a concrete 
example, consider again the three images in Figure 1. The /3-model mass estimates yield 
X = l.l7, 1.02 and 0.97 from top to bottom whereas the scaling law results are ^ = 1.05, 
1.00 and 0.99, respectively. 

The "noise" added by the measured values of (3 has several sources. Recent dynamics play 
a role. Clusters are not in exact hydrostatic balance, particularly in their outer parts. The 
present density and temperature structure can be perturbed by prior mergers and accretion. 
Geometry also plays a role. The cluster gas distribution is, in general, ellipsoidal rather 
than spherical and the measured values of (3 are obtained from projected, two-dimensional 
images of the three-dimensional density distribution. Finally, values of (3 derived from 
surface brightness fits are sensitive to contamination from foreground/background sources, 
choice of cluster center, as well as image quality and technical aspects of data reduction 
procedures. These concerns are compounded when one considers the basic fact that (3 is a 
measure of the derivative of the logarithm of the surface brightness. 

The scaling law method is superior to the /3-model because of its smaller variance. Its 
accuracy is also remarkably insensitive to the dynamical state of clusters and the cosmo- 
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logical background. Its main drawback is the reliance on numerical experiments to provide 
the normalization, rio(S c ), of equation (8) which, in turn, depends on the particular struc- 
ture formation model under scrutiny, and may depend as well on the numerical method 
used in the investigation. Regarding the latter, we are encouraged by the good agreement 
shown between two independent codes used in this study. We anticipate future studies 
by groups employing independent numerical methods will address the robustness of the 
normalization ri (5 c ). We also find room for optimism in the rather modest sensitivity to 
fi displayed by the normalization and scatter in mass estimates shown in Table 5. This 
insensitivity may be rooted in the fact that the scaling laws merely reflect the condition of 
virial equilibrium within 5 C ~(9(10 3 ). Overall, these results show that it is possible to use 
X-ray spectroscopy to estimate cluster masses with an rms accuracy substantially better 
than 20%. 

6. Discussion and Conclusions. 

The results described in the previous sections agree well with those of Schindler (1995), 
who found biases and variance of a few tens of percent or less in the "normal cluster" sam- 
ple derived from the numerical simulations of Schindler & Bohringer (1993) and Schindler 
& Miiller (1993). The simulations and analysis methods in those works differ in many 
respects from those used here. In particular, the Eulerian gas dynamics scheme adopted in 
their study captures shocks more accurately than the SPH technique used here, while the 
adaptive nature of the SPH smoothing kernel provides better spatial and mass resolution 
in high density regions. Given the significant differences between the independent simu- 
lation algorithms used in these studies, we find the degree of qualitative and quantitative 
agreement rather encouraging. 

Tsai, Katz & Bertschinger (1994) also examined the accuracy of the /9-model mass es- 
timates applied to an SPH simulation which included radiative cooling for the gas, but 
which excluded the effects of galaxy formation and feedback. At 1 Mpc from their cluster 
center, where 5 C ~ 700, they found that the /3-model overestimated the mass by ~ 25%, 
typical of the uncertainties found in our analysis. They also found large underestimates 
at radii near the core radii of their surface brightness fits (~200 kpc) which arose, in part, 
from the presence of a strong gradient in the gas temperature near the center. However, 
at these high density contrasts, their results are compromised by the artificially strong 
concentration of baryonic material near the center. 

The improved variance of the scaling law estimates over the /9-model is a mixed blessing. 
A serious concern is the sensitivity of the slope and normalization of the rs c — Tx relation 
to the assumed cluster physics and the underlying cosmological model. We find very 
modest sensitivity to fio near 5 C ~ 10 3 (Table 5) and find that the ejection models do not 
differ substantially from the infall sample. We suspect, and this suspicion is supported 
by collisionless simulations of cluster formation (Crone et al. 1994), that the impact of 
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changing the initial perturbation spectrum will be less than that of varying O , but this 
issue remains to be investigated in detail. 

What should then an observer with X-ray data do to estimate a cluster's mass? It depends 
on whether it is more important to minimize the bias or the variance in the mass estimate. 
For example, when comparing a set of X-ray derived masses to values derived using an 
independent method (e.g., weak gravitational lensing), minimizing the bias in the X-ray 
binding masses would perhaps be most important, and the /3-model approach preferred. 
On the other hand, if one were looking for the slope of correlations between an observable 
cluster property (e.g., optical luminosity or velocity dispersion) and binding mass, then 
minimizing the variance would be more important, and the scaling law method preferred. 

To summarize, we have used an ensemble of 58 gas dynamics cluster simulations to investi- 
gate the accuracy of binding mass estimates based on the hydrostatic, isothermal /9-model 
and on the temperature-mass relation. We have also analysed the velocity and tempera- 
ture fields of the numerical models to address the questions of hydrostatic equilibrium and 
isothermality. A summary of our main results follows. 

• Within a radius, rsoo, where the cluster mean interior density is 500p c , the gas velocity field 
is extremely quiet (Table 3 and Figure 3), validating the basic assumption of hydrostatic 
support by thermal pressure. Despite local variations in temperature due to ongoing merger 
events (Figure 1), the radially averaged gas temperature is nearly isothermal within rsoo 
in the = 1 sample. The O = 0.2 clusters exhibit a moderate, negative radial temperature 
gradient. 

• The standard /3-model mass estimator (eq. 4) is nearly unbiased and has a modest scatter 
in regions where the mean estimated density contrast is in the range 500 < 5^, st < 2500. 
For example, at 8 e c st = 1000, the mean value of M est /M true is 0.94 (1.15) with standard 
deviation 0.23 (0.19) for the = 1 (O = 0.2) sample. The bias and scatter both increase 
with cluster radius (decreasing 8 C ). The bias increases because the true density profiles 
are steeper than the assumed isothermal value, while the dispersion increases because of 
the longer dynamical timescales characteristic of larger radii. 

• We find a strong correlation between rs c , the radius encompassing a mean density contrast 
S c , and Tx, the broad beam, emission- weighted gas temperature. This scaling — a reflec- 
tion of the similarity between clusters of different mass and their near virial equilibrium 
state — can be used to generate mass estimates with smaller variance than that of the 
/3-model. The degree of scatter is surprisingly small. At 5 C = 1000, the standard deviation 
is 0.11 (0.12) for the = 1 (Oo = 0.2) sample. The larger dispersion in the (3- model method 
arises because two parameters, f3 and T x , contribute independent sources of error whereas 
the scaling law method incurs error from only one parameter, Tx- 

The results from the experiments reported here and those from other experiments cited 
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above practically rule out the possibility of large systematic errors in the mass determina- 
tion of galaxy clusters. The large baryon fractions measured in clusters therefore remain 
difficult to reconcile with standard primordial nucleosynthesis in an O = 1 universe. 
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Table 1 

Summary of Model Notation 



Label 


N 


Description 


Reference 


EJ 


14 


= 1, as = 0.59, wind ejection included 


Metzler & Evrard (1994,5) 


2F 


14 


same as EJ without winds 




NFW 


6 


= 1, a 8 = 0.63, Tree/SPH code 


Navarro et al. (1995a) 


EdS 


8 


= 1, as = 0.59, no feedback 


Evrard et al. (1993) 


F12 


8 


O = 0.2, A o = 0.8, as = 1, same IC's as EdS 




0p2 


8 


O = 0.2, A o = 0, a 8 = 1, 





Table 2 

Synthetic Observation Parameters 



Parameter value 



source redshift z 0.04 
exposure time t exp (sec) 7200 

background noise level <S/v a 3 x 10 -4 

background subtraction level Si, a 3 x 10~ 4 



minimum surface brightness in fit S m i n a 3 x 10 
a units : ROSAT cnts s -1 arcmin -2 . 

Table 3 

Mean Mach Numbers within rsoo 



Sample 


(v r )/c s 


(vl)/cl 


EJ 


0.001 ±0.016 


0.041 ±0.010 


2F 


-0.022 ±0.022 


0.069 ±0.018 


NFW 


-0.080 ±0.119 


0.261 ±0.077 


Eds 


-0.008 ±0.033 


0.112 ±0.019 


F12 


-0.012 ±0.004 


0.022 ±0.009 


Op2 


-0.005 ±0.014 


0.045 ±0.020 
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Table 4 

/?-Model Accuracy 





= 1 




O = 0.2 




5 est 


x a 


°x 


X a 


°x 


100 


1.46 


0.53 


2.08 


0.50 






n if, 

u.ou 




u.oo 


500 


1.02 


0.29 


1.34 


0.24 


1000 


0.94 


0.23 


1.15 


0.19 


2500 


0.87 


0.16 


1.00 


0.14 


X = M 


est j^jtrue 












Table 5 








Scaling Law Accuracy 






= 1 




O = 0.2 




5 C 


rio(5 c ) a 


ox 


rio(£ c ) a 


a x 


100 


4.89 


0.20 


4.78 


0.20 


250 


3.37 


0.18 


3.31 


0.16 


500 


2.48 


0.15 


2.48 


0.14 


1000 


1.79 


0.11 


1.87 


0.12 


2500 


1.11 


0.08 


1.25 


0.10 



Defined by equation (8), units are Mpc (/i = 0.5). 
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Figure Captions 



Figure 1. Synthetic ROSAT X-ray surface brightness (left column), emission weighted 
temperature (center column) and emission weighted, projected velocity field (right column) 
for three orthogonal projections of a single cluster from the EdS sample. The cluster is at an 
assumed redshift 2 = 0.04 and bars in the figure show a 10' angular scale. In the grey scale 
contours, the dark (or light) bands are logarithmically spaced by factors of 10° 6 in surface 
brightness, 10 2 in temperature, with the third from minimum dark band representing 
10~ 3 cnts s _1 arcmin -2 and 10 7 K, respectively. The velocity vectors are spaced every 2', 
and scaled such that l' = 1000 km s -1 . 

Figure 2. Projected azimuthally averaged X-ray surface brightness (top) and emission 
weighted temperature (bottom) for the cluster shown in Figure 1. The top to bottom lines 
correspond to the top to bottom projections, with data for the middle and bottom projec- 
tions displaced by factors of 10 (2) and 100 (4) for the surface brightness (temperature) 
with respect to the top projection. 

Figure 3. Radial Mach numbers (v r )/c s (left) and (r 2 )/c 2 (right) derived from the gas 
sound speed c 2 (r) =5kT(r) /3fj,m p and the first and second moments of the radial velocity 
field for all the runs in the ensemble. All quantities are local values measured in radial 
shells. The overdensity 5 C is used as a radial coordinate; note the inverted axis. The 
dashed line shows 5 C = 500, our conservative estimate for the boundary of the hydrostatic 
region. 

Figure 4. Scaling between cluster size, as measured by rsoo, and emission weighted 
temperature for all the models. Symbol types correspond to different models, as shown in 
the legend. The data are well fit by equation (5). Each model appears three times, from 
three orthogonal projections. 

Figure 5. Three dimensional temperature profiles for all the clusters in the ensemble. 
The temperature in radial bins is expressed in terms of the average, emission weighted 
temperature Tx and radius is normalized to the cluster size rsoo- 

Figure 6. Accuracy of the /3-model mass estimates as a function of the estimated den- 
sity contrast 8^ st , equation (6), for the ensemble. Each model appears three times, from 
orthogonal projections. 

Figure 7. Histograms of the estimated-to-true mass ratios derived from the /3-model 
evaluated at three different estimated density contrasts. Dashed vertical lines show an 
error of ±40%. 

Figure 8. X-ray surface brightness maps of the six worst underestimates from the = 1 
ensemble. Values of the estimated-to-true mass ratio are shown above each panel. Within 
each panel, the light and bold circles represent the true and estimated values of rsoo, 
respectively. Strongly bimodal or complex images usually result in poor /3-model mass 
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estimates. 



Figure 9. Histogram of the estimated-to-true mass ratios from the /3-model for the = 1 
combined sample at S^, st = bOO. 

Figure 10. Histogram of the estimated-to-true mass ratios from the scaling law method 
for the = 1 combined sample at 5 C = 500. The distribution is unbiased by construction. 

Figure 11. Scatter plot of the deviations in (3 and Tx (defined in the text) for all the 
models in the ensemble. Point styles are the same as those used in Figure 4. 
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